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Neutrinoless hadron Lepton Number Violating (LNV) decays can be induced by virtual 

Majorana neutrino, "which in turn indubitably show the Majorana nature of neutrinos. Many 

three-body LNV processes and Lepton Flavour Violating (LFV) processes have been studied 

extensively in theory and experiment. In this current paper, as a supplement, we study 

"^ ■ 75 four-body LNV (LFV) processes from heavy pseudoscalar B and D decays. Most of 

these processes have not been studied in theory and searched for in experiment, while they 

Q^i may have sizable decay rates. Because the four-body decay modes have the same vertexes 

^5 I and mixing parameters with three-body cases, so their branching fraction are comparable 

with the corresponding three-body decays. We calculate their decay widths and branching 

fractions with current bounds on heavy Majorana neutrino mixing parameters and hoped 

\N . these four-body decay processes could be probed in experiment. 
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00 ! 1 Introduction 
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f^ . In the Standard Model (SM), neutrinos are strictly massless, but non-zero neutrino masses have 

m ■ 

been detected in experiment [1-4]. So considering the physics of neutrinos, extension of the SM 
is surely needed. But until now, the nature of neutrinos is still puzzled, we do not know whether 
they are Dirac or Majorana particles. So before we determine the way how to extend the physics 
of SM, we have to clear the Dirac or Majorana type of neutrinos. 

There is strong theoretical motivation for Majorana mass term to exist since it could naturally 
explain the smallness of the observed neutrino masses [5, 6]. We all know, though not coming 
from first principle, the SM conserves the lepton number, but Majorana mass term violates 
lepton number by two units (AL = 2), thus the neutrinoless hadron LNV decays with like sign 
dilepton final state are crucial for the existence of Majorana neutrinos. The possible Lepton 
Flavor Violating (LFV) meson decays could be induced either by Majorana neutrino or neutrino 
oscillation in which case the neutrino is a Dirac neutrino. However, neutrino oscillation at loop 
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level would be suppressed by powers of -^^ and thus could not bring the branching fraction to 
an observation level. As a result any direct observation of LFV (LNV) process indicates the 
existence of Majorana neutrino. 

There were many attempts to determine the Majorana nature of neutrinos by studying the 
LNV and LFV processes. Since the neutrinos in the final state are undetectable by the detectors, 
therefore neutrinoless processes are preferred, e.g., the neutrinoless double /3 nuclei decay (Oz//3/3) 
has long been advocated as a premier demonstration of possible Majorana nature of neutrinos 
[7, 8]; the Majorana neutrino exchange r lepton three-body or four-body decays [9-11]; the LNV 
process pp — )■ £^£^ + X oi pp ^ £^£^jj at LHC [12, 13]; the top-quark or W-boson four-body 
decay [14, 15]; the LNV or LFV meson decays with like sign dilepton in the final state [16-21], 
et al. 

Recently, Atre et al. [22] have studied K, D, Dg and B decays via a fourth massive Majorana 
neutrino, they found if the exchanged Majorana neutrino is resonant, that means the exchanged 
Majorana neutrino is on mass-shell, then the corresponding branching fractions can be tremen- 
dously enhanced by several orders, which may be reached by the current experiments. Inspired 
by the effect of resonant neutrino, various AL = 2 meson three body decays M^ — )■ ^^I'^^^Mg" 
have been studied in lecture [22-27] and in lecture [28] the authors discussed four-body decays 
B -^ DUtt. 

In the experiment, some of these LNV (LFV) processes have been searched for. For example, 
Fermilab E791 Collaboration reported their results of searching for the LNV and LFV decays of 
D^ into 3 and 4-bodies, they presented upper limits on the branching fractions at 90% confidence 
level (CL) [29]. Recently, using 772 x 10^ BB pairs accumulated at T(4S') resonance with the 
same CL, the Belle Collaboration set the upper limits on the LNV (LFV) B^ — t- D^i^i + decays 
[30]. Using a sample of 471 it 3 million BB events, the BABAR Collaboration searched for the 
LNV processes B — t- K^(TT^)£'^i^ and placed upper limits on their branching fractions also 
with 90% CL [31]. The LHCb CoUaboration, using 0.41 fb'^ of data collected with the LHCb 
detector in proton-proton collisions at a center-of-mass energy of 7 TeV, reported their upper 
limits on the branching fractions of B~ decays to D^*'^ p^ p^ , ir^p^p^ and D^p^p^ at 95% 
CL. They also searched for the 4-body decay B — )■ D^n^p^ p^ and set upper limit for the first 
time on its branching fraction [32]. The experimental situation of searching for the LNV and 
LFV processes can be found in the papers Refs. [11, 33]. Though these LNV and LFV processes 
are still non-observation, but the upper limits for branching fraction have been obtained, which 
also in turn limit the mixing parameters between Majorana neutrino and charged lepton. 



Though lots of LNV (LFV) processes have been studied in experiment or in theory, but 
there are stiU many channels which have not been considered. Especially the four-body LNV 
(LFV) meson decays, most of them are still missing in literature, and some of this kind channels 
may have sizable branching fractions and may be reachable in current experiment. LNV four- 
body decays can also offer complementary information about the masses and heavy mixings 
of such a heavy (resonant) Majorana neutrino, so they need to be studied. In this paper, we 
study 75 four-body LNV (LFV) processes of dilepton decays B[D) — )■ Miiii2M2, where Mi 
stand for a pseudoscalar meson, M2 can be a pseudoscalar or a vector meson, ^1(^2) = e, /U. 
These AL = 2 LNV (LFV) 4-body meson decays are induced by a Majorana neutrino, and 
the possible lowest order diagrams are illustrated in Figure 1 (a-b). Some processes, like the 
decays of B^ — )■ D^-f^^^Mg", where Mg" stand for vr^, p^, K^ , D^ or D^, only have the 
Feynmann diagram shown in Figure 1 (a); but some decays, like B^ — )■ D^£^i2M ~, where 
M ^ denote D^ or D^ , have both decay modes shown in Figure 1 (a) and (b). In Figure 1 (a), 
if the Majorana neutrino mass lies between the range from a few hundred MeV up to 4.4 GeV 
(since it is heavy, it may be a fourth generation neutrino), then the neutrino mass could be on 
mass-shell (resonance), and the corresponding decay rate is much enhanced due to the effect of 
neutrino-resonance, so its contribution will be much larger than the one of neutrino-exchange 
diagram in Figure 1 (b), which is suitable for a continuous neutrino mass. So in this paper, 
we only consider the contribution of neutrino-resonance diagram Figure 1 (a) and ignore the 
contribution of neutrino-exchange diagram Figure 1 (b). We also eliminate contribution of the 
interference between the two diagrams. 

There are two key points to calculate these 4-body decay modes. One is the choice of the 
mixing parameters, since most of these Majorana neutrino induced 4-body decay modes do not 
have experiment results, we cannot extract the mixing parameters through these decays. So 
we followed the method by Atre et al, they determine the parameters by experimental data 
[22]. We choose the strongest constrains which were got from the current data as input in our 
paper [22, 26], purpose is to ensure exactly. The other is the calculation of the hadronic matrix 
element between initial meson B{D) and final meson Mi, we use the Mandelstam formalism [34] 
where the hadronic matrix element is described as a overlapping integral over the wave functions 
of the initial and final states [35]. The wave functions are obtained by solving the relativistic 
Bethe-Salpeter (BS) equation [36]. 

This paper is organized as follows, in section 2, we outline the formulas of the transition 
matrix element. In section 3, we give the details of how to calculate the hadronic matrix 
element. In section 4, we give the results and draw the branching fraction of heavy meson 



4-body decays as a function of the heavy neutrino mass. 



2 Theoretical Details 



The leading order Feynman diagrams for the LNV (LFV) 4-body decays of heavy meson M: 



M{P)^Mi (Pi )it{P2)4 {P^)M^ {Pa) 



il) 



are shown in Figure 1 (a-b). Here M is the pseudoscalar B or D with momentum P, two charged 
leptons £'^ , £2 have momentum P2 and P3, pseudoscalar meson Mi with momentum Pi denotes 
vr, K or D, meson M2 with momentum P4 can be a pseudoscalar meson vr, K, D and DgCt al, 
or a vector meson p, K* , et al. 

Such LNV (LFV) process can occur through a Majorana neutrino, and the vertex between 
this Majorana neutrino and charged lepton is beyond the SM. Following previous studies [22, 37], 
assuming there is only one heavy Majorana neutrino which may be a fourth generation neutrino. 
It can be kinematically accessible in the range we are interested in. Then the gauge interaction 
lagrangian responsible for the LNV (LFV) decay can be written as: 



C = -^W+ V V;,Ni^^^PL£ + h.c. 



V2 



(2) 



where Pl = 2(1 ~ Ts)) ^4 is the mass eigenstate of the fourth generation Majorana neutrino, 
V£4 is the mixing matrix between the charged lepton i and heavy Majorana neutrino A^4. 
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Figure 1: Feynman diagram of the four-body decay of heavy meson 
The transition amplitude for the 4-body decay M(P)-^Mi(Pi)^|(P2)£^(P3)M2"(P4) shown 
in Figure 1 (a) can be written as: 



M 



9%^<i.V, 



<ii<i2 "qsqi 
8M^ 



(Afi(Pi)|gi7^(l - 75)g2|M(P)) X M^^u X (M2(P4)|^37"(1 " 75)«4|0), (3) 



where we have delete the momentum dependence in the propagator of W boson since it is 
much smaller than the W mass; g is the weak coupling constant; Vq^q^ (Vq.^q^) is the Cabibbo- 
Kobayashi-Maskawa (CKM) matrix element between quarks qi and 52 (qs and (74); Ai^u is the 
transition amplitude of the leptonic part. 

As mentioned in Introduction, we only consider the contribution of the diagram in Fig- 
ure 1 (a), where the Majorana neutrino is on mass shell, and the effective narrow- width approxi- 
mation relating the resonant contribution can enhance the decay rate substantially. In this case, 
according to Ref. [22, 26], the leptonic matrix element Aifj,u is given as: 



52 

M^u = yV£i4V;?24"1-4 



UlJ^'yuPR'^2 Ui-f^-f^PRU2 



ql. -m\ + iTN^niA q'^ -ml + iTNirn4 



(4) 



where Vg4 is the mixing parameter between the heavy Majorana neutrino and charged lepton; 
Pr = f (1 + 75); Qn4 is ^^^ momentum of heavy Majorana neutrino (g^^ is the case of exchange 
the two final charged leptons); 771,4 is the mass of the heavy Majorana neutrino; r7V4 is the total 
decay width of the heavy neutrino. 

Mesons M and Mi are pseudoscalar mesons, the corresponding hadronic matrix element in 
Eq. (3) can be described as a function of form factors: 

(Mi(Pi)|^i7^(l - 75)g2|M(P)) = P^iU + /-) + P^{f+ - /-), (5) 

where the method to calculate the form factors /_|-, /_ will be given in next section. 

The last part (M2|/i2|0) in Eq. (3) is related to the decay constant of the meson M2. If M2 
is a pseudoscalar with momentum P4, we have the following relation: 

(M2(P4)|g37"(l - 75)94|0) = iFM.P:^, (6) 

where Fm2 is decay constant of meson M2. If M2 is a vector with momentum P4 and polarization 
vector e, then the corresponding relation becomes to: 

(M2(P4,e)|g37'(l - 75)^7410) = M2FM,e\ (7) 

here we use the same symbol M2 to denote the meson and its mass. 

Combining Eq. (4), Eq. (5) and Eq. (6), we rewrite the decay amplitude Eq. (3) in the case 
of meson M2 as a pseudoscalar: 

M = 2GpVe^4Ve^4Vq-^q^Vq^q^FM2m4, 

^- \ FMf+ + /-) + A A(/+ - /-) , nnf+ + n + AA(/+ - n ] p ... 

where Gp is Fermi constant. If meson M2 is a vector, we just replace ^^4 with M2 /in numerator 
in Eq. (8). With the numerical values of form factors /_[_ and /_ obtained in next section, the 
calculation of this decay amplitude is straightforward. 



3 Hadronic transition matrix element 



In order to calculate the hadronic matrix element and get the numerical value of form factors 
/+, /_, we choose the Mandelstam formalism [34], where the transition amplitude between two 
mesons is described as a overlapping integral over the Bethe-Salpeter wave functions of initial 
and final mesons [35]. Using this method but with further instantaneous approximation [38], in 
center of mass system of initial meson, in leading order, we write the hadronic matrix element 
as [39]: 



(Mi(Pi)lgi7^(l-75)g2|M(P)) 



dq 



■Tr 



V^Pi^ (9*1)7^(1 - 75)^V(q) 



M 



(9) 



(27r)3 

where P and Pi are the momenta of initial and final mesons; M in denominator is the mass 
of initial meson; q is relative momentum between quark and antiquark inside the initial meson; 



Qi 



q + 



m2 



—r-, -f IS the relative momentum inside the final meson Mi, m\ (rriL) is mass of 

antiquark (quark) in final meson Mi, r is three dimension momentum of meson Mi; (p^^ is the 



positive wave function for a meson in the BS method; for the final state, we have use the symbol 

lP^J 70- 



'PpI 



loifV 



Table 1: Mass of quark in unit of GeV. 



quark 


b 


c 


s 


d 


u 


mass(GeV) 


4.96 


1.62 


0.5 


0.311 


0.305 



In the BS method, the positive wave function 99++ for a pseudoscalar meson can be written 
as [40]: 



^P 



A[B + j^ + 4^C + 



.r 



M 



^75, 



(10) 



where q±_ = (0, g), and 



A = 


M 
2 


, ,^ , , ,^mi+m2 

h{q) + f2{q) ■ 

UJi + UJ2 


B = 


Wl + UJ2 


C = 


mi — m2 


■miL02 + m2UJi ' 


D = 


Wl + ijJ2 



[ID 



miUJ2 + m2UJi 

In Eq. (11), mi and m,2 are the masses of quark and antiquark inside the meson, and we list 
their values in Table 1; Ui is defined as uJi = \/m,f + (p, i = 1,2; fi{q) and /2(^ are the wave 
function of the meson. 
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With Eq. (10) and Eq. (11), we take the integral on the right side of Eq. (9), then the form 
factor /-I-, /_ can be expressed as: 

1 fTi T2 M - El, 



■'^ 2 Im ^ Ml ^ 



M 



-T, 



/- 



1 /Ti T2 M + Ei 

2 [m Ml M 



T3 



(12) 



where Mi and Ei = \/ M\ + r^ are the mass and energy of final meson Mi ; and 



Ti 
T2 
T3 

h 



(27r)3 
d?q 

1 f (fq 



AAiAti, 
4AiAt2, 
^4AiAt3\q\cose, 



f\ J (27r) 

mi2 Cqi-q BDi / ^ ^ , "212 ^ 



mil +mi2 
-CCi if + 



Ml Ml 



mil + mi2 



mi2 _^2 



"111 +"1-12 



77112 \ mio E] 

-qi-q]-BBi-DDi ^ T^^i • 9, 



-1 



"1-12 



"^11 +"1-12 



mil +"1-12 
-CiMi 



inn +mi2 Ml 



^^^ -BDiEi - DDif, 



mil +mi2 



-Ci-BiD-BDi^-^ + ^[2qi.q + 



mi2 _2 



(13) 



Ml Ml Ml V '^ ' mii + mi2 

where Ai, Bi, Ci and Di have the same definitions with Eq. (11) but replace the parameters 
by the one of final pseudoscalar. 

Numerical values of wave functions fi{q) and f2iQ) can be obtained by solving the coupled 
Salpeter equations [40]: 



(M - 2a;i) 



(M + 2a;i) 



fM + f2{q) 



7T1-1 
^1 



/ 



dk 1 



(2vr: 



i3,,,2 



fM - f2{q) 



mi 

wi 



+f2ik)miiui 
dk 1 



(27r)3 ojf 



iVs-V,)[fiik)ml 
-{Vs + V,)fi{k)k-q} 
{Vs-V,)\fi{k)ml 



/2(fc)mi6^i] - {Vs + V,)fiik)k ■ q} . (14) 



Where we have chosen the Cornell potential, which is a linear potential plus a single gluon 
exchange reduced vector potential, and in momentum space the expression is: 

1 



Vs{q) 
Vv{q) 



- + Vo) 6^{q} + — 7^—^72 : 

a J TT^ (g^ + a^)^ 

; as{q) 



37r2 (g2 + Q,2) 
127r 1 

^log(a+T# 



(15) 



A^ • 



where a = e = 2.71828; A = 0.21 GeV^ is the string constant; a = 0.06 GeV is a parameter to 
avoid the infrared divergence; the QCD scale Aqcd = 0.27 GeV characterize the running strong 
couphng constant a^; the constant Vq is a parameter by hand in potential model to match the 
experimental dada, whose values for different mesons are listed in Table. 2. 

Table 2: Parameters Vn in unit of GeV 



meson 


B 


D 


K 


vr 


Vo 


-0.091 


-0.375 


-0.962 


-0.999 



With these parameters, we solved the full Salpeter equation Eq. (14), and obtained the 
numerical values of wave functions fi{q) and f2{<i) for pseudoscalar mesons B, D, K and vr, at 
the same time the meson masses of these pseudoscalar mesons are also obtained which consist 
with experimental data. 

4 Numerical Results and Discussions 

Besides the parameters appearing in potential, there are other parameters whose values needed 
to be assigned. We choose the CKM matrix elements [41]: Vud = 0.974, Vus = 0.225, Vcd = 0.230, 
Vcs = 0.973, Vcb = 40.6 x 10^'^, Vub = 3.89 x 10~'^. The decay constants of pseudoscalar and 
vector mesons used in our calculation are listed in Table 3. 

Table 3: Decay constants Fm2 of pseudoscalar and vector mesons in unit of MeV. 



meson 


vr 


P 


K 


K* 


D 


Ds 


Fm, 


130.4 [41] 


220 [42] 


156.1 [41] 


217 [42] 


222.6 [43] 


260 [41] 



The key step to calculate the decay widths and branching fractions of LNV (LFV) heavy me- 
son decays is to determine the limits on the mixing parameters |V£j4Vg24l ^^d the heavy neutrino 
mass 7724 in Eq. (4). Following the approaches of Refs. [22, 26], we take the mixing parameter V^4 
and the mass 777.4 as phenomenological parameters. Since the mixing parameters are common 
factors, we choose some decay modes have the same {Vij^^Vi^^l under our consideration and have 
mixing parameters numerical upper bounds in experiment, thus we extract the numerical values 
of mixing parameters from these processes, details can be found in Refs. [22, 26]. We choose 
the strongest constrains on mixing as input in this paper, purpose is to ensure exactly. For 
the value of 7774, since we only consider the case when the heavy neutrino is on mass shell, so 
we determine the mass of neutrino by kinematics. With numerical values of mixing parameter 



and neutrino mass 1114, the neutrino total decay width r^r^ is calculated, which is summing all 
possible decay channels of Majorana neutrino at the mass 7714 [22]. So in our calculation, T]\f_^ is 
mass and mixing parameter dependent, not a fixed value. 

With these parameters and the limits on mixing parameters, 75 LNV (LFV) decay widths 
and branching fractions of the heavy mesons D^, D^, B^, and B^ are calculated. Among 
these processes, there are some channels where the meson Mi is a light meson, vr or K. We 
should point that since we have made instantaneous approach to Bethe-Salpeter equation, so 
the calculation of the hadronic matrix element including a light meson may not be accurate as 
the heavy meson case. But since all these decays are beyond the SM, precise calculation is not 
the issue, so we also give the results including these decays. In the calculation of decay rate, we 
perform a Monte Carlo sampling of the branching fractions and the mass of heavy neutrino, i.e., 
we calculate the excluded region of the branching fractions as a function of the heavy neutrino 
mass m^, and we plot the results in Figures 2-7. Where the regions inside and above the curve 
are excluded by current experiment data, where below the curve is allowed in theory. 

The curve is not smooth, this phenomenon is caused by two reasons. First, we choose 
different mixing parameters {Vij^^Vi^il due to different ranges of heavy neutrino mass 777.4. Since 
the current limits on mixing parameters are related to heavy Majorana neutrino mass, that is, 
according to different neutrino mass range, we choose different LNV (LFV) processes to get the 
strongest constrains on mixing parameters. For example, in process B^ ^ D^ e^ e'^ M2 , we choose 
three processes K^ — t- e"^e+7r^, D^ — )■ e+e+vr^ and B^ — t- e+e+vr^ to limit |144p. Second, as 
announced above, the value of neutrino total decay width T^^ is mass and mixing parameter 
dependent, whose value is changed along with the changes of different mixing parameter |Vg4| and 
neutrino mass 7714. As a result of the piecewise mixing parameters, the branching fractions are 
also piecewise as a function of the neutrino mass. This different value choice of mixing parameters 
may be the main reason to cause the differences between our results from the one in Ref. [28, 44], 
they also calculated the branching fractions of B^ — )■ D^e^e^ir^ and B^ — )■ D^fi^fi^Tr^. 

There is another phenomenon seems unusual in the results of some branching fractions. For 
example in Fig. 5 (b), we show the branching fraction for the decay mode D^ — t- K^e^fi^K^. 
At two edges of the curve, which are the locations of the allowed smallest and largest neutrino 
masses, the values of the branching fractions are very small. This behavior of small rates is not 
unusual actually, it happens due to the restriction of the phase space, the very small kinematic 
phase space at edges lead to these small branching fractions. 

Because some 4-body decays of mesons have broader phase space than the corresponding 



3-body processes and the resonance neutrino mass is determined kinematically, so one advantage 
of these 4-body decays is that we can detect much wider range of Majorana neutrino mass. For 
example, we can study the heavy neutrino if its mass is in the range of 2 GeV — )■ 4 GeV through 
the 3-body decay B^ — )■ e~e^D^ [26], while through the 4-body decay B^ — )■ L'"e^e+7r~ or 
B^ — 7- D^e^e^TT^ , we can reach the range of possible neutrino mass from 0.2 GeV to 3.4 GeV. 
Another advantage of these 4-body decay is that the branching fraction is not small comparing 
with the corresponding 3-body decay [22, 26], because they are in the same order, in some cases, 
they have same vertexes, that is, they have same mixing parameters |V£^4Vg24l ^^d CKM matrix 
elements. 

We have mentioned that the dominant determinant of the branching fractions comes from 
the mixing parameter \Vi^iVi,^4\, which are limited by the current experimental data. Besides 
these parameters, there are another important parameters, CKM matrix elements, which are 
also determinant factors to the values of branching fractions. We note that if the final mesons 
are D~ and vr", there are two decay modes, B^ — )■ D~if£2'^~ and B^ — )■ vr^^j'^^^D^. In 
the first decay mode, the CKM matrix elements are |T4bKtd|^) while for the second one are 
jKfoV'cdl^, as {VubVcdl"^ /iVcbVudl"^ ~ 4 X 10""^, so we ignore the decay B^ — ^ TT^ifi^^^ ^.nd its 
interference with B^ — )■ D^^^^^^vr^. For the same reason we only consider the contribution of 
decay D^ — )■ K^ifi2'^^ and ignore the decay mode D^ — )■ vr^^j'^^^i^^. 

In B decay modes, from the figures we can see that, if the heavy neutrino mass is located 
in the range of 0.8 — 3 GeV, the maximal branching fraction of -B — t- DliH.2'^ would approach 
10~®. This outcome is very impressive because the correspond experimental result is Br(B^ — )■ 
D^fi^H^TT^) < 1.5 X 10^® [32]. For the cases of D decays, we showed in the figures that if the 
heavy neutrino mass is between 0.5 GeV and 1.5 GeV, the branching fraction of D — t- Kiii2T^ 
almost reach 10^^. This number is bigger than B decay modes, and just below the correspond 
experimental upper limit Br(D^ — )■ K^fi^fi^-K^) < 3.9 x 10~^ [29]. Finally, we point out 
that, in calculations of the decay modes B^ — )■ 7r'^£^^^M2', we lack the information of mixing 
parameter |T44V^4| when neutrino mass 1714 > 4 GeV, so the results in Figure 6 (b) are given by 
set {VeiV^^l = in these cases, that is, there are no predictions when neutrino mass rn.4 is larger 
than 4 GeV in Figure 6 (b). 

In some particular channels, there is additional contribution coming from intermediate 
mesons resonance [44]. For instance, in the decay channel i?+— T-D^/i^/i+vr^, except the CKM 
favored diagram drawn in Figure 1 (a), there is another CKM dis- favored diagram (see Fig- 
ure 1 (b) in Ref. [44]), where the two final mesons can be induced by a intermediate resonance 
15*^(2010), in range of 2.1 GeV < 1714 < 3.3 GeV, the intermediate resonance D*~(2010) may 
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results in a sizable contribution in decay B^ ^ D^ fi^ fi^ tt , but we do not consider these cases. 

In conclusion, we extended the existing literature to consider the 4-body LNV (LFV) rare 
decays of heavy mesons B and D, since the 4-body decays share the same vertexes and mixing 
parameters as well as the CKM matrix elements with the corresponding 3-body decays, relative 
large branching fractions which are comparable with the 3-body decays are obtained, which can 
be searched for in the existing and forth coming experiments. 
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